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Final  Report 

Grant#:  FA9550-09- 1-0403 
David  Bortz 

Department  of  Applied  Mathematics 
University  of  Colorado,  Boulder 

Abstract  We  proposed  a  novel  approach  which  employs  random  sampling  to  generate  an  ac¬ 
curate  non-uniform  mesh  for  numerically  solving  Partial  Differential  Equation  Boundary  Value 
Problems  (PDE-BVP’s).  From  a  uniform  probability  distribution  U  over  a  ID  domain,  we  con¬ 
sidered  M  discretizations  of  size  N  where  M>#.  The  statistical  moments  of  the  solutions  to  a 
given  BVP  on  each  of  the  M  ultra-sparse  meshes  provide  insight  into  identifying  highly  accurate 
non-uniform  meshes.  We  used  the  pointwise  mean  and  variance  of  the  coarse-grid  solutions  to 
construct  a  mapping  Q{x)  from  uniformly  to  non-uniformly  spaced  mesh-points.  The  error  con¬ 
vergence  properties  of  the  approximate  solution  to  the  PDE-BVP  on  the  non-uniform  mesh  are 
superior  to  a  uniform  mesh  for  a  certain  class  of  BVP’s.  In  particular,  the  method  works  well 
for  BVP’s  with  locally  non-smooth  solutions.  We  fully  developed  a  framework  for  studying  the 
sampled  sparse-mesh  solutions  and  provided  numerical  evidence  for  the  utility  of  this  approach  as 
applied  to  a  set  of  example  BVP’s. 

Summary  Over  the  duration  of  this  grant,  while  developing  our  SMRT  methodology  for  solving 
BVP-PDEs,  the  core  of  our  research  efforts  have  include  the  following:  substantial  refinement 
to  our  algorithm,  extension  of  the  algorithm  to  higher  dimensions,  and  establishing  the  theoretical 
well-posedness  of  our  approach  [3,4].  All  of  these  topics  are  linked  by  a  desire  to  efficiently  exploit 
the  high  paralellizability  of  our  approach  and  future  implementation  on  massively  parallel  multi¬ 
core  technologies.  Lastly,  we  have  also  been  invited  to  contribute  a  review  article  on  computing 
on  GPU’s  to  SIAM  Review  [5].  The  focus  of  this  effort  is  one  type  of  computation  which  is 
substantially  accelerated  on  GPU’s. 

We  now  give  a  brief  summary  of  our  progress. 

Scandalously  Paralellizable  Mesh  Generation  The  PI  and  his  collaborator  are  developing  an 
SMRT  framework  to  generate  non-uniform  meshes  for  solving  PDE’s  [3,4,5].  These  discretizations 
can  offer  superior  solution  accuracy  and  convergence  properties  to  that  of  uniform  spacing.  We 
offer  a  brief  overview  of  our  proposed  algorithm  as  well  as  the  establishment  of  a  preliminary 
theoretical  framework  [3].  Also,  in  [4]  we  extended  resutls  in  [2]  to  the  identification  of  (9  using 
an  optimization  technique  using  results  from  probability  theory.  However,  we  discovered  that  the 
approximation  technique  described  below  was  substantially  more  efficient. 

We  consider  a  monotonically  non-decreasing  function  (9 :  I  — ^  I  which  is  absolutely  continuous 
on  a  finite  number  of  compact  subsets  of  I  and  restricted  at  the  endpoints  to  (9(0)  =  0,  (9(1)  =  1. 
The  purpose  of  the  function  (9  is  to  map  the  uniformly  spaced  mesh  to  a  non-uniformly  spaced 
one.  The  goal  is  to  develop  a  strategy  for  identifying  a  (9  such  that,  e.g.,  the  approximate  solution 
to  the  Poisson  problem 

u{Q(x))=f{Q(x))  s.t.  u(Q(0))=A;  u(Q(l))=B, 
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has  convergence  properties  (in  n)  superior  to  a  uniform  spacing.The  core  of  our  approach  is  to 
identify  Q  via  a  sparse  stochastic  approximation.  We  repeatedly  sample  from  a  distribution  P  and 
then  use  pointwise  statistical  moments  of  the  coarse  solutions  to  generate  the  desired  non-uniform 
mesh  function  Q.  Naturally,  different  classes  of  problems  call  for  different  strategies  for  generating 
Q.  Our  results,  however,  suggest  that  a  more  generalizable  strategy  may  exist.  Before  presenting 
our  conclusions,  we  briefly  establish  some  notation. 

Let  p  be  a  function  taking  a  point  q  e  I  and  a  random  vector  of  length  n,  and  mapping  them  to 
a  single  random  variable 

p&X(n)(P))  EJs:[{£/(X(B)(P))}Jf=JXw=^]  .  (1) 

The  function  U  takes  a  discretization  of  the  domain  and  solves  the  BVP.  The  operator  denotes 
expectation  with  respect  to  a  uniform  distribution  on  { 1 , . . . ,  n)  where  the  distribution  of  the  index 
random  variable  K  and  { -  }  ^  denotes  the  kth  element  of  a  vector.  We  note  that  this  function  returns 
a  random  variable  for  each  L,.  Let  the  pointwise  mean  of  p  be  defined  for  £,  e  I  as 

Ml)  =Ep  [Ejt  [{tf(Xw(P))}s._JX(i)  =5]]  .  (2) 

The  pointwise  variance  of  p  is  defined  for  £  e  I  as 

v©  =  Vp  [Ejr  [{E/(Xw(P))}s._t|Z(i)  =5]]  ,  (3) 

where  Vp  denotes  variance  with  respect  to  P.  Ek  denotes  expectation  with  respect  to  U  { 1 . . . . ,  11 } , 
the  distribution  of  the  index  random  variable  K,  and  {■}/<  denotes  the  kth  element  of  a  vector. 

Answers  to  the  critical  questions  for  this  approach  are  depicted  below 

For  each  candidate  Q,  how  many  sample  sparse  grids  need  to  be  generated?  The  relation¬ 
ship  between  the  mesh  size  n  and  the  number  of  samples  m  is  non-trivial.  and  Figure  1  illustrates 
this  by  depicting  the  error  in  v  (relative  to  v  computed  with  m  =  3000  sampled  from  a  uniform 
distribution  on  I)  for  a  range  of  n  and  m  values.  For  a  given  n,  though,  we  do  note  that  the  error 
in  the  v  computation  is  decreasing.  In  Figure  2  we  depict  the  number  of  samples  of  vector  size  n 
which  are  needed  to  ensure  three  digits  of  accuracy  in  estimating  the  variance.  Since  the  number 
was  consistently  below  1000  over  a  range  of  n,  we  let  m  =  15000  in  all  subsequent  simulations 
(unless  otherwise  specified). 

In  what  way  do  the  random  solutions  converge  to  the  actual  solution?  For  a  conventional 
finite  difference  discretization,  we  would  consider  the  error  E  in  the  solution 

\\E(Q,  *!)||  =  \\u(Q(x°n))-U(Q(x°n))\\ 

=  |Ae(x0)  -/Q(xo))  | 

-  A0(x0)  X0(x°)  » 
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v(U( I))  error  vs.  m  and  n 


Figure  1:  log10  of  the  error  in  the  computation  of  v  (sampling  from  a  uniform  distribution  on  I)  as 
a  function  of  m  and  n.  Note  the  general  downward  trend  along  both  the  m  and  n  axes. 
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Figure  2:  For  each  n,  the  vertical  axis  reflects  the  number  of  samples  needed  to  compute  the 
variance  with  3  digits  of  accuracy  relative  to  v  (sampling  from  uniform  distribution  on  I)  with 
m  =  3000. 
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which  is  bounded  above  by  the  spectral  radius  of  the  inverse  of  the  finite  difference  operator  A 


and  a  truncation  error  t 


G(i 


GO'S) 

For  the  non-uniform  three-point-stencil  approximating  the  second 
derivative,  the  truncation  error  is  0(max/{  \hk\).  For  our  development,  we  consider  a  probabilistic 
version  of  this  error,  with  the  following  conditions. 

CONDITION  Cl .  For  a  given  P,  the  spectrum  of  A ^  ^  is  bounded  in  [0, 1]. 

CONDITION  C2.  For  a  given  P,  the  truncation  error  induced  by  a  finite  difference  approxima¬ 
tion  to  the  second  derivative  is  first  order  in  the  largest  step- size  h. 

Theorem  1 .  Under  Condition  Cl  and  C2,  the  expected  error  converges  pointwise  to  zero. 

See  [3]  for  support  of  these  conditions  as  well  as  a  proof  of  the  theorem. 


How  should  Q  be  constructed?  The  function  Q  is  created  using  the  statistical  moments  of 
the  sampled  sparse-mesh  solutions  and  based  on  results  in  [1],  For  the  problems  with  second 
derivatives  we  define  Q  as 


Q(x) 


(*)> 


where 

91 M  =  f  yV(i;;f/(X(„)(P))|<«; . 

and  the  superscript  —  1  is  an  inverse  function  operator.  Essentially,  this  definition  will  pile  up  points 
in  regions  with  a  steep  solution  in  an  effort  to  provide  higher  order  accuracy  for  the  nonuniform 
second  derivative  discretization. 

For  the  problem  with  a  second  power  of  the  first  derivative,  we  define  Q  as 


Q(x) 


'<12  (')' 
.92(1). 


(*)> 


where 

Jo 

and  v  is  defined  above.  Evidence  for  improvement  in  error  convergence  is  depicted  in  Figures  3-4 


We  hypothesize  that  the  reason  c/\  (jc)  works  well  is  that  the  p'  may  converge  faster  than  p.  We 
also  hypothesize  that  the  function  <72  (t)  works  well  because  the  second  derivative  (when  cast  as 
the  local  curvature)  is  inversely  proportional  to  the  local  variance  of  a  random  variable  (a  result 
which  is  well  known  in  the  semi-parametric  nonlinear  regression  literature).  Essentially,  while  the 
p"  may  not  converge  quickly,  the  product  p"v  does.  We  also  found  that  multiplication  by  an  extra  v 
dramatically  improves  the  computed  Q ,  though  an  explanation  is  not  immediately  clear.  A  deeper 
understanding  of  the  spectrum  of  AX(;i)(/>)  and  how  it  depends  upon  the  choice  of  P  will  be  essential 
to  explaining  the  efficiency  of  <72 (*)•  We  plan  to  explore  both  of  these  issues  in  a  future  paper  [4]. 


Acknowledgment/Disclaimer  This  work  was  sponsored  (in  part)  by  the  Air  Force  Office  of 
Scientific  Research,  USAF,  under  grant/contract  number  FA9550-09- 1-0403.  The  views  and  con¬ 
clusions  contained  herein  are  those  of  the  authors  and  should  not  be  interpreted  as  necessarily 
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Error  Convergence  For  Each  Mapping 


Mesh  Size  N 

Figure  3:  Error  convergence  for  uniformly  and  non-uniformly  spaced  points  for  the  steady-state 
Hamilton- Jacobi  BYP. 


Figure  4:  Error  convergence  of  the  different  mesh  mappings  for  the  singular  BVP. 
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